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ABSTRACT 


In^Zhis  •  paper  we  describees  procedure  for  solving  the  inverse  scattering 

problem  in  ultrasonic  imaging  .^Althcugji'the  derivation  of  the  method  is  based 

M'/i t- 1 

on  the  Helmholtz,  equation  model. 

A  ^ 


A  +  k2(l+f)u 


it  is  applicable  to  any  model  for  vdiich  the  spacial  sound  pressure  u  *  u(r,k) 
satisfies  an  asynptotic  equality  of  the  form 

(2)  *(k)  s  ^  log  u(r,k)  |*d  -  J  F(r>ds  +  0(k'°)  ,  k  -  « 

rs 

where  k  is  proportional  to  the  frequency  P  demotes  the  ray  path  along  which 
the  pressure  wave  travels  from  the  source  point  rg  to  the  detector  point  r^ 
and  o  is  a  positive  constant.  The  method  is  based  on  predicting  ♦(*)  ■  /p  F(r)ds 
via  a  rational  function  procedure ,  using  several  values  $(kj)  ,  $0^)  ...  , 

The  method  is  illustrated  for  the  case  of  the  Rytov  approximation  to  (1) ,  in  vdiich 
case  the  paths  P  are  straight  lines.  A  perturbation  method  for  correcting  for 
curved  ray  paths  is  also  described.  The  algorithm  can  also  be  modified  to  image 
materials  with  more  complicated  frequency  dependent  attenuation.  Examples  of  images 
reconstructed  from  computer  simulated  data  with  and  without  Gaussian  additive 
noise  are  given.  The  beneficial  effect  of  a  noise  tolerant  first  norm  data  fitting 
algorithm  in  improving  image  quality  is  shorn. 


Key  Words:  Ultrasound,  tomography,  imaging,  Helmholtz  equation,  rational  function, 
attenuation,  multifrequency,  diffraction,  refraction. 
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The  methods  derived  in  this  paper  are  based  on  the  Helirholtz  differential 
equation  node! 


(1.1) 


V^u  +  k2(l  +  f)u 


where  u  denotes  the  spacial  sound  pressure,  v  is  the  Laplacian  operator  in 
3  dimensions ,  and 


(1.2) 


f© 
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In  Eq.  (1.2),  r  -  (x.y.z)  ,  Cq  is  the  speed  of  soisid  in  the  region  (usually  a 
liquid)  surrounding  the  body  B  ,  k  *  u/Cq  ,  where  u  m  2tt  x  (frequency)  and  c(r) 
is  the  "speed”  of  sound  at  the  point  r  in  B  .  We  assune  that  B  is  a  subset 
of  a  volume  V  ,  and  that  sotnd  sources  and  detectors  are  located  on  8V  ,  the 
boundary  of  V  .  We  have  displayed  the  word  "speed"  in  quotes  above,  since  in  the 
derivation  of  Eq.  (1.1)  c(r)  is  indeed  the  speed  of  sound  at  the  point  r  ,  whereas , 
in  real  life,  materials  have  attenuation,  and  therefore  2rr  f  >  C  .  In  this  paper 
we  derive 
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a  numerical  al-.-  ‘  **n  for  reconstructing  the  function 


a.3> 


F  -  /"ITT 


in  V  ,  based  on  the  model  (1.1). 

The  fact  that  the  function  c(r)  turns  out  to  be  ccrplex,  as  well  as  other 
criticisms  of  (1.1)  have  raised  questions  with  regards  to  the  validity  of  the 
model  (1.1),  and  this  has  led  several  authors  [1,4]  as  well  as  sane  of  us  [3] 
to  derive  other  models.  However,  under  either  plane  wave  or  point  source  exitaticn , 
all  of  these  newly  derived  models  including  (1.1)  have  the  caiman  feature  that 


(1.4) 


1 

IK 


log  u^(r) 


F  ds  +  o(l)  ,  k  -*•  ® 
L 


where  ^(r)  is  the  spacial  sound  pressure  and  L  is  a  path  in  V  along  tfttLch 
the  pressure  wave  travels  from  its  source  point  rg  to  the  detector  point  ?d  . 

While  other  inversion  techniques  are  based  on  the  particular  model,  the  method  of 
this  paper  applies  to  any  such  model. 

Although  f  and  are  not  C®  functions  of  r  in  practice,  and  they  are 
therefore  difficult  to  confute,  we  show  that  the  function 

(1.5)  $00  -  log 

is  a  very  smooth  and  slowly  varying  function  of  k  man  interval  [k^,*)  of  the 
real  line,  where  kQ  >  0  .  Indeed,  on  tk0,»3  ,  the  function  «>00  satisfies  the  xelatia 


(1.6)  $(k)  -  J  F  ds  +  0(k'6)  ,  k  -  •  , 

where  6  is  a  positive  constant  depending  only  an  f  .  In  addition,  although  we 


2 


cannot  give  a  couplet*  proof,  we  expect  that  +(k)  is  analytic  and  bounded  in  8 
sector 


(1.7)  S(kg,0)  -  {k  €  tt:  |k|  £  kg  ,  |arg  k|  <  0} 


and  where  kg  and  6  are  positive  constants .  Recently  discovered  results  by 
approximation  theory  [8]  have  shown  that  such  a  function  can  be  very  accurately 
approximated  on  [kg,®]  by  a  rational  function  of  k  .  Thus  if  $(k)  is  known 
only  on  a  subinterval  of  [kg,®]  ,  we  can  expect  the  "p-algprithn'Csee  e.g.  [12]) 
to  yield  a  very  accurate  approximation  to  $(®)  ®  f^Fds. 

To  this  end,  we  have  simulated  sane  tests  for  the  case  of  the  Rytov  approxi¬ 
mation  to  the  solution  of  (1.1),  taking 


(1.8) 


f(r) 


e 


and  when  u£  is  a  plane  wave.  In  this  case  it  is  possible  [6,7]  to  explicitly 
evaluate  the  Rytov  approximation,  Ryt[u^(r)]  and  the  paths  L  are  straight 
lines  L  parallel  to  the  incoming  field  u£  .  By  computing 


(1.9) 


4^00 


Ryt .  (u^(r) } 


at  9  equi -spaced  points  in  the  range  of  frequencies  between  1  megahertz  and 
4  megahertz,  the  p  algorithn  enables  us  to  accurately  compute 


This  single  algorithn  works  well*  if  the  mmbers  $(k)  can  be  ccoputed 
accurate  to  4  or  more  significant  figures,  but  breaks  ckxn  at  3..  ’In  practice,  we  expect 
that  we  will  be  able  to  measure  ♦(k)  accurate  to  3  significant  figures.  To  this  end,  we 
have  introduced  an  "i  -adjustment"  to  our  algorithm,  which  in  effect  computes  the 
limit  of  several  rational  function  interpolations  of  $(k)  ,  and  then  evaluates 
the  minimum  of  these  limits.  This  has  the  effect  of  under-erphasizing  the 
larger  errors  of  these  limits,  and  if  the  $(k)  are  accurate  to  3  significant  figures, 
we  are  able  to  compute  $(»)  accurate  to  3  significant  figures. 

In  order  to  reconstruct  an  image  fr cm  the  #(») ,  we  use  standard  X-ray  CT 
algorithns  [2J  which  assume  that  the  paths  L  are  straight  lines  L  .  While  this 
approximation  is  usually  not  unduly  pessimistic  for  the  case  of  sound  waves  in  tissue, 

It  Is  not  exact.  To  this  end,  we  have  introduced  a  numerical  scheme  for  correcting  the 
T  vfriich  we  constructed  under  the  assumption  that  the  paths  are  straight  lines.  This 
method  1*  based  an  a  discretization  of  the  equations  recently  derived  in  [7] .  While 
we  have  included  this  algorithm  for  ccnpleteness,  we  have  not  vet  tested  it  in  practice. 

In  the  final  section  of  this  paper  we  illustrate  a  reconstruction  via  the 
p-algariths  of  the  function 

(1.11)  f<r)  -  l  «'b^'rj>2 

3-1 

using  simulated  data  obtained  via  the  Rytcrv  approximation.  We  also  Illustrate  the 
result  when  0.1%  noise  was  added  to  the  Rytov  data,  and  finally,  we  illustrate  the 
result  when  our  ^-adjustment  algorithn  Is  applied  to  the  p-algorithm  results 
computed  using  the  noisy  data.  It  is  possible  that  we  could  have  obtained  better  re¬ 
sults  by  first  applying  the  smoothing  algorithm  to  the  data,  and  then  applying 

the  rational  extrapolation  algorithm.  We  shall  investigate  this  procedure  in  the  future. 
- ~  — ~ — — 

In  our  simulation  we  are  able  to  distinguish  t*ro  objects  of  the  form  (1.8) 
separated  by  2  wavelengths. 
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In  our  model  (1.1)  we  have  assumed  Chet  f  is  independent  of  m  •  We 
could  still  recover  f  if  its  functional  form  were  faotci  explicitly,  such  as,  if 


(l.U)  fCfc.r)  -  o(5)k6  +  B(r>  . 

by  suitable  normalization  of  lqg  u^(r)  ,  e.g. ,  by  altering  (1.5)  to 

aia  *00  -  log  . 

This,  then  would  enable  us  to  recover 

(1.13)  *(-)  -  |  ds  . 

Next,  we  could  recover  S(r)  by  sorpllng 

(1.14)  ¥(k)  -  k6[*(k)  -  ♦(•)] 

-  [  *±g©d*  ,  k*-  . 

4  2a*(?) 
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2.  INTEGRAL  EQUATION  POftflJLATIQN  OF  THE  HELMHOLTZ  EQUATION 


By  means  of  Green's  theorem,  one  can  derive  the  integral  equation 


(2.1)  ^(r)  -  u°<?)  +k2  |j|  (^(r-P)f(P)u(r’)dV(r*)  . 


which  is  equivalent  to  (1.1),  where 


(2.2) 


"TJttt 


and  where  u£(r)  denotes  the  input  special  pressure  wave.  In  practice  one  usually 
has  point  sources  of  the  form 


(2.3) 


»|r-r.| 
e _ _ 

Att  lr-r_  I 


in  which  r  denotes  a  source  point,  air  hough  the  plane  wave  form 
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(2.4) 


is  sometimes  also  vised,  where 


(2.5) 


<vyv  •  w  - k  *  £ 
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ASYMPTOTIC  ESTIMATE  OF  log 

Let  be  defined  by  either  Eqs.  (2.3)  or  (2.4).  The  theory  of  geometric 
optics  then  predicts  that  (see  e.g.  [3  p.  134)) 

(3.1)  ^-3-  -  «p{ik[f  ,f+7ds  +  o(l)]>  ,  k»«  . 

Jl 

where  in  the  case  of  u£  defined  by  (2 . 3) ,  L  denotes  the  ray  path  emanating 
fran  rg  to  the  point  r^  ,  vhile  in  the  case  of  defined  by  (2.4),  1 

starts  out  in  the  direction  of  k  at  some  point  an  the  boundary  of  V  and  ends 
up  at  rd  . 

Let  us  now  assure  that  f  €  Lip^(V)  ,  for  some  a  >  0  .  We  believe  that  in 
applications  we  must  always  have  o  >  1/2  ,  although  for  purposes  of  deriving  our 
algorithm,  any  positive  value  of  o  will  suffice.  Then,  upon  substituting  (3.1), 
with  ?d  replaced  by  r  into  (2.1)  and  assuming  L  to  be  any  C1  path,  vr  find 
that  the  o(l)  term  in  (3.1)  can  be  replaced  by  0(k“°^2)  .  That  is, 

(3.2)  3^-  -  «p(lk[f  .T+T  ds  +  0(k*°^2)])  .  k-»  . 

\<*s>  ‘ 

This  result  may  be  established  via  a  procedure  similar  to  that  applied  to  the 
equations  (2.2)  and  (1.7)  of  110],  with  o  >  0  ,  n  >  0  and  m  +  n  »  a  . 
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4.  ANALYTIC  PROPERTIES  OF  TOE  SOLUTION 

In  this  section  vie  study  the  analytic  properties  of  the  solution  of 

(2.1)  as  a  function  of  k  ,  in  order  to  show  that  we  can  perform  accurate 
approximation  of  the  function  $  defined  in  (1.5)  on  [kQ,»]  . 

We  shall  shew  that  with  u£  defined  by  (2.3),  the  solution  u^Cr)  of 

(2.1)  is  a  ratio  of  two  entire  functions  of  k  .  The  poles  of  this  function 
occur  at  the  points  k  ■  k^  for  which  the  homogeneous  equation  corresponding 
to  (2.1)  has  solutions.  We  then  shew  that  these  ntrbers  cannot  be  real. 

3 

Let  V  be  a  bounded,  region  in  IT  ,  and  let  f  €  Lip^  (V)  ,  such  that 

3 

Jm  f  >  0  a.e.  on  V  ,  and  such  that  f  =  0  on  It  -  V  .  Hence 

(4.1)  +  k2[l+f(r)]Uk(r)  -  0  ,  r  €  V 

(4.2)  ^(r)  +  k\(x)  -  0  ,  r  6B3  -  V  . 

3 

Now  Wilcox  [11]  has  shewn  that  if  Im  k  >  0  on  It  •  V  ,  then  any  solution 
u  of  (4.2)  must  satisfy 

2 

(4.3)  |  |  |i  -  iku|  dR(r)  -  0  ,  R  -  - 

r*R 

and 

(4.4)  J  |u)2  dS(r)  -*•  c(0  <  c  <  ■>)  ,  R  -  *  . 

r-R 

*y  Green's  theorer,  we  have 

(4.5)  |  uv^udVCr)  j  |?uj2dV(r)  »  J  u|pdS(r) 

r<R  r<R  r-R 
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that  is, 

(4.6)  -k2  |  (1+f)  |uJ2dV(r)  +  |  |7u|2dV(r)  -  |  u^dS(r)  . 

•  r<R  r<R  u-R 

Let  us  now  assume  that  k  is  real.  Then,  by  taking  Imaginary  parts  of  (4.6), 
we  get 


(4.7) 


(Im  f)  |u|2dV(r)  »  Im  J  u  |H  dS(r) 

r*R 

-  Im  J  u[||  -  iku]dS(r)  +  Im(-ik)  J  |u|2dS(r)  . 
r*R  r»R 


By  applying  Schwartz's  Inequality  to  the  first  Integral  on  the  right  hand 
side,  letting  R  ®  ,  using  (4.3)  and  (4.4),  and  recalling  our  assumption 

that  Im  f  >  0  a.e.  on  V  ,  we  get 


(4.8)  0  <  k2  |  (Im  f)  |u|2dV(r)  -  -  kc  <  0  . 

V 

This  contradiction  proves 


UEMMA  4.1:  If  Im  f  >  0  a.e.  on  V  ,  then  the  equation 


(4.9) 


?  r  -ik|r-r' I 

u^(r)  -  Yf-  j  5 - ><*<?> 

v  4n  j  r-r '  j 


cannot  have  my  nontrivial  solutions  for  k  real. 


Indeed,  attempts  at  extending  Lenina  4.1  to  canplex  values  of  k  suggest 
that  if  we  set 
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* 


|  f(r)dV(r)  -  peie 
V 

|V|  -  f  dV(r)  , 
Jv 


where 


(A.  11)  p  >  0  ,  0  <  e  <  ir  ,  p/|V|  <  1/2 

and  if  6*  in  the  range  -ir/8  <  6'  <  0  is  defined  by 

(4.12)  tan  26’  -  -  ■  -P~siP_e.— 

f  v  |+p  cos  e 

then  the  equation  (4.9)  has  no  nontrivial  solutions  if  k  *  |k|eia  satisfies 

(4.13)  -  e’  <  a  <  *  -  6‘  . 

Although  we  cannot  prove  this,  we  expect  that  the  following  assumption  is  probably 
valid  if  the  equation  (4.1)  has  no  solutions  for  which  the  ray  paths  are 
trapped  in  V: 

ASSUMPTION  4.2:  If  Im  f  >  0  in  V  and  if  (4.10)  and  (4.11)  are  satisfied 
then  there  exists  a  sector 

(4.14)  S(k0,6)  -  {ke  C:  |k|  >kQ  ,  |argk|  <6} 


where  and  6  are  positive  numbers,  such  that  whenever  (4.9)  has  a  nontrivial 
solution  u^  then  k  4  SCk^.e)  . 


•j 


LBfft  4.3:  7!-  '  'lution  u  of  Eq.  (2.1)  may  be  represented  In  the  farm 

(4.15)  “|c®  * 


where  far  each  fixed  r  ,  X  and  Y  are  entire  functions  of  k  . 

PROOF;  We  carry  out  the  proof  only  for  the  case  of  u£  defined  by  (2-3). 
The  proof  in  the  case  of  (2.4)  is  similar,  and  we  emit  it. 

Instead  of  (2.1),  we  consider  the  equation 


u(k,a,r) 


ialr-r 


4tt  r-r 


(4.16) 


+  k 


ia|r-r 


ATT  )r-rs } 


8  f(P)u(k,a,r’)dV(r’) 


We  know  frem  the  theory  of  Fredholm  integral  equations  (see  e.g.  [4])  that  for 

every  fixed  a  €  C  and  r  •  (x.y.z)  ,  the  solution  u  of  Eq.  (4.16)  is  a  ratio 

2 

of  two  entire  functions  of  k  ,  and  hence  of  k  ,  i.e., 


(4.17) 


Y(k,a,r) 


where  X  and  Y  in  (4.17)  are  entire  functions  of  k  .  Inspection  of  the 

proof  of  this  result  for  the  case  of  (4.12)  shows  that  for  fixed  k  and  r  , 
X  and  Y  in  (4.17)  are  also  entire  functions  of  a  .  Hence  for  fixed 
r  ,  X(k,k,r)  and  Y(k,k,r)  are  entire  fvnctiais  of  k  . 


11 


IQtft  4.4:  If  f  e  Lip0(V)  ,  if  u£  is  defined  by  either  Eqs.  (2.3)  or 
(2.4),  and  if  Assumption  4.2  holds,  then  the  function  $  defined  by  (1.5) 
is  analytic  and  bounded  in  the  sector  S(kg,6)  ,  where  it  satisfies  the  relation 


(4.18) 


*00 


ds  4  0(k‘o/2) 


k  •*  •>  in  SO^.e)  . 


PROOF:  By  Assumption  4.2  there  are  no  nontrivial  solutions  of  Eq.  (4.9), 

with  k  in  S(kg,6)  .  Hence  the  function  Y(k,r)  in  Eq.  (4.15)  does  not 
vanish  if  k  e  SOfyS)  .  Hence  *(k)  is  analytic  in  SOfye)  .  Eq.  (4.18) 
therefore  follows  from  Eq.  (3.2). 


5.  RATIONAL  APPROXIMATION  PREDICTION 


In  this*  section  we  shall  show  that  we  can  "predict" 

(5.1)  $(*)  ■  j  »4  +  f  ds 

using  the  Thiele  algorithm  for  constructing  a  rational  function  which  interpolates 
$(k)  .  We  first  prove 


LBfrlA  5.1:  Let  $(k)  be  defined  as  in  Lemma  4.2.  Corresponding  to  every  integer 
N  >  0  ,  there  exists  a  rational  function 


(5.2) 


r„00 


c/  +  c/'1  ♦ 
kM  +  a/'1  + 


such  that 


(5.3) 


I $00  -  rN(k)  |  <  C  e 


W2 


kp  <  k  ^  ® 


where  c  and  C  are  positive  constants  independent  of  N  and  tdiere  kg  is  de¬ 
fined  as  in  Laima  4.1.  Indeed,  for  any  c  >  0  ,  we  may  take 


(5.4) 


c 


,1100% 

(~r> 


1/2 


c 


*#iere  o  is  defined  as  in  (3.2)  and  6  as  in  (4.10).  Then  C  •  C(e)  . 


PROOF:  In  [9]  a  rational  function  of  z  similar  to  (5.2)  was  constructed,  for 
approximating  a  function  $(z)  on  -1  <  *  <  1  ,  where  $  was  analytic  on  the  unit 
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disc.  A  simp!: 


#aticn  of  the  argment  In  [9]  yield*  (5.3)  and  (5.4). 


The  Inequality  (5.3)  tells  us  that  rational  functions  db  a  considerably  better 
job  of  approximating  $(k)  than  polynomials  In  1/k  .  Fbr  example,  if 
I*00 -♦(»)!  -  ckm°f2  ,  c  +  0  ,  then  no  natter  hew  the  constants  *>Q,\ . b^,  are  selected. 


(5.5) 


max 

*0^- 


H 

1*00  -  l  b, 

J-0  J 


where  C  is  a  positive  constant  vhich  depends  neither  on  Ihe  nor  on  N  . 

In  practice,  the  constant  o  Is  not  known  a  priori.  Nevertheless  by  the 
linear  methods  of  [8,9],  (5.3)  is  still  satisfied  for  same  positive  constants  c 
and  C  . 

Vie  note  also,  from  (5.1),  (5.2)  and  (5.3),  that 


(5.6) 


!♦(»>  -  *N(»)  I 


This  fact,  and  the  fact  that  Lenina  5.1  tells  us  that  ♦(k)  is  'Very  nearly”  a 
rational  function,  suggests  that  the  Thiele  rational  function  algorithm  should 
provide  an  effective  procedure  for  determining  Cq  . 

Given  the  data 

2N 

(5.7)  {k.  ,  *(k,)} 

J  J  j-0 

where  the  kj  are  distinct  and  real,  the  Thiele  Algorithm  for  constructing  an 
rN(k)  of  tl*  form  of  Eq.  (5.2)  is  described  by  the  following  equations: 
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6.  1  IMPROVEMENT  OF  NOISY  DATA 

The  algorithm  (5.8)  for  predicting  $(“)  has  been  tested  on  the  function  (1.9), 
under  the  assurption  that  we  can  compute  $00  in  the  range  of  frequencies  between 
1  and  A  megahertz.  Indeed,  taking  N  as  low  as  3  and  evaluating  $(k) 
accurate  to  7  decimals  we  can  compute  the  integral  in  (1.10)  to  6  decimals 
accuracy.  Our  test  runs  have  shewn  however,  that  **iile  (5.8)  still  yields  reliable 
results  when  the  nunbers  $(kj)  are  rounded  to  4  decimals  accuracy,  it  fails 
when  they  are  rounded  to  3  .  We  have  foind  that  this  problen  can  be  remedied  by 
taking  the  (usual)  minimun  of  these.  That  is  if 


<cT> 


m«  1 


are  M  such  predictions  obtained  by  applying  (5.8)  to  M  different  sets  of  data 
of  the  form  (5.7),  the  i1  minimum  of  the  Cq  is  the  number  a  which 


minimizes 


(6.1) 


I  itf4-  .i 

nri. 


3h  this  way,  we  are  able  to  use  noisy  data  of  the  form  (5.7)  vhich  is  accurate  to 
only  3  decimals  to  predict  $(*)  ■  a  accurate  to  3  decimals. 

In  the  case  when  the  c£j  are  real,  it  is  easily  seen  that  minftS (a)  •  min^SCcQ11^) 
Hence,  even  for  general  complex  ,  it  is  convenient  to  replace  the  exact  solu¬ 

tion  a  of  (6.1)  by  the  approximation,  a  *  o+i6  ,  there  a  and  B  are  real,  and 


vhere 


(6.2) 


m  He 


ej*0:  S(1V)  =  l  I**  4’)-«'liSa)f<  C<a'h 


»  -  *=  S<»<n  -=  I  |>C0W.»'|,S®(>a  . 
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7.  CORRECTION  FOR  RAY  REFRACTION 

Our  experience  has  shown  that  the  ray  paths  L  of  soind  waves  in  tissue  are  usually 
very  nearly  straight  lines  L  .  We  may  thus  use  the  equation  derived  in  [7]  to 
derive  a  correction  on  the  f  *  F  ,  where  F  is  constructed  under  the  assumption 
that  the  paths  L  are  straight  lines  L  .  To  this  end,  if  r '  •  (0,0,0)  , 
rd  -  (*.0,0) 

(7.1)  Tl  »  |  F  ds  •  F(x,e4>(x),ciKx))  A  +  e V (x)*  +  eV (x)^  dx  , 

that  is.  L  is  the  path 

(7.2)  L  -  {(x,y,z):  y  -  e$(x)  ,  z  *  etKx)  ,  0  <  x  <  *} 

and 

(7.3)  F  -  vTTF 

In  the  notation  of  [7] ,  we  write 

(7.  A)  F  -  1  +  ch  , 


i.e. 


(7.5) 


eh 


»T+T  -  1  . 


In  ultrasonic  tonography  |c(r)  -  c0l  /  |cQ|  <  1/10  ,  so  that  |f|  <  1/5  ; 

hence  assuming  that  |h|  <  1  ,  we  may  as  sine  the  relation  |e|  <  1/10  . 

In  order  to  reconstruct  F  via  straight  line  methods,  we  would  like  to  know 


(7.6) 


To  this  end,  the  following  equations  are  derived  in  17] 

Sm 

(7.7)  e  •  J  e2  ^  \  $’2  +  j  V2)  dx  +  0(e3)  , 


where 


(7.8) 


hy(*.y.*) |y  .  i  .  0 
h2(x,y,z)  |  m  z  m  0 


as  well  as  the  equations 

(7.9)  ♦"  -  hy  ,  -  \  . 

which  are  solved,  subject  to  the  boundary  conditions 

(7.10)  <K0)  -  ♦(*)  -  ♦(O)  -  *(*>  -  0  . 

ye  now  apply  these  results  to  our  case. 

Assisting  that  the  paths  are  straight  lines,  we  can  reconstruct  an  F  ,  such 

that 
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(7.11) 


fl  . 

\  •  Lro*t  ; 


Me  would  like  to  knew  Fq  ,  such  that 


(7.12) 


Then,  we  have 


-  F(x,0,0) 


'  t'o*  : 


-  F(x,0,0) 


(7.13) 


Fi-e 


tdiere,  upon  ignoring  terms  of  order  t  ,  and  taking  e  ■  1  (implying  a  corresponding 
reduction  in  the  magnitudes  of  the  corresponding  quantities) , 


(7.14) 


F0  +  7  [*,(x)2  +  *’(x)2J  , 


where  by  (7.8).  (7.9)  and  (7.10) 


(7.15) 


*’(x)  -  if  (t-i)F  (t,0,0)dt  +  [X  F(t,0,0)dt 

4  Jo  y  '0  y 

ip'(x)  *  |  jZ  (t-£)Fz(t.0,0)dt  +  j*  Fl(t.0.0)dt 


Let  us  new  derive  a  numerical  scheme  for  making  the  correction  (7.14).  To 
this  end  we  assume  that  F  is  fooun  on  a  3-dimensional  pixel  array  with  side  h  . 


whose  centers  are  located  at  the  points 


I  J  K  A  U 

(7.16)  {(x4,yvz.)}  -  (((i-i)h  ,  a4)h  .  (k-i)h}A 

^  3  ^  i-1  j-1  k-1  2  2  4  1.1  iml  ^ 


I  J  K 


We  then  use  the  approximations 
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(7.17) 


r-y  =  s  _  l.J-M.K  1.1-1 

Fijk  -  3y  «L*rV  2k  — 


»  _  Fij.k+1  *  Fi1.k-1 

^jk  —  2k  — 


f  ry(t,yj .V*  *7hi?jk  +  hjJ%k  +  7h fy k  i  Aljk 


C  ?  (t'yJ  •  V*  *  7  h  ^jk  +  *2  +  7  h  %k  5  B: 


i-1 

I 

m*2 


sijk 


£x  I 

f  j0  (t*V  By  ft.yj'V*  *  £j,  {<m‘7)h-V^k5  c; 


x  npl 


mjk  ’  jk 


£ 

^|0X  <**V  yt.yyV*  *  til  {0"'7)h'V^k5  v 


■JO  get 


7.18) 


?ijk  *  Fijk  *  7  ^Aijk  *  Cjk)2  +  ^ijk  *  Djk)2j 


After  making  such  an  "x-sweep"  for  each  j  and  k  ,  we  then  make  a  similar 
’y-sweep",  and  repeat,  witil  the  changes  in  the  numbers  F,.v  become  negligible. 


8.  AN  EXAMPLE 


In  this  Motion  we  illustrate  the  application  of  the  formulas  developed  in 
Sections  5*6  of  the  paper  via  a  simulation  carried  out  for  the  case  of  the  Rytov 
model.  The  Rytov  approximation  u(r)  *  Ryt.{u(r)}  is  often  made  for  the  solution 
to  (1.1),  where 


(8.1) 


Ryt.{u(r)} 


iE*r+tl 

e  * 


and  where 


(8.2)  U^(r) 


k2  f  f  [  expdklr*?'  |  +  Ht»(r*r,))f (?*)<¥(?*) 

|r-?| 

IT 


Our  simulation  was  carried  out  for  the  case  of 


--  .2 


(8.3) 


f(?)  - 


4  -blr-r»| 

I  *  J 

j-1 


since  it  is  possible  to  explicitly  evaluate  W  when  f  has  the  form  (8.3)  (10). 
In  (8.3)  the  points  Fj  are  located  at  •  (O.yj  ,0)  ,  there  in  wits  of  milli- 
msters,  y!  *  3.75,  y2  *  ®  •  ?3  *  *2*5  »  7 4  •  -4.375  .  The  timber  b  was  choew 
so  that  each  Gaussian 


-bfr-rJ2 

a  J 


had  a  full  width  at  half  maxima  equal  to  In,  l.a. ,  about  2  wavelengths  at 
a  frequency  of  2.5  W  (megahertz),  and  souid  speed  Cq  ■  1500  m/sec. 

All  reconstructions  are  based  on  the  formula 


t 


-  |rg  a 

• 

where  L  Is  i  straight  line  path  [10],  and  where  the  numbers  (8.4)  were  obtained  by 
using  the  algorithms  (5. 8) -(5. 10)  by  taking  9  values  of  k  corresponding  to 
9  equi-spaced  frequencies  in  the  range  of  1  to  4  MR  . 

Figure  8.1  is  a  picture  of  the  exact  object  digitized  to  a  41  x  41  pixel 
array,  each  pixel  having  dimension  .3125  tun  x  .3125  ton. 

The  pictures  in  Figures  8.2,  8.3  and  8.4  (as  well  as  Figure  8.1)  were  obtained 
from  reconstructions  carried  out  in  the  Ultrasound  laboratory  at  the  University  of 
Utah,  and  based  on  the  methods  of  this  paper.  The  pixel  array  in  each  dimension 
is  the  same  as  in  Figure  8.1.  The  projection  data  was  obtained  from  30  views  over 
180  degrees.  The  line  of  detection  was  positioned  at  6  cm  from  the  center  of 
rotation. 

The  picture  in  Figure  8.2  was  obtained  by  applying  the  algorithms  (5-7) -(5. 10) 
to  data  {k. ,  2W.  (r)/(ik.)}^  accurate  to  7  decimals.  It  deviates  from  the 

J  J  j-1  _ 

exact  object  at  each  end,  since  there  the  value  of  f(r)  in  (8.3)  was  very  small, 
so  that  the  resulting  significant  figure  accuracy  in  the  data  was  less  than  4 
places.  The  rational  fixiction  algorithm  was  therefore  less  accurate  at  these  ends. 

The  picture  in  Figure  6.3  is  a  reconstruction  obtained  as  in  Figure  8.2,  but 
by  adding  random  noise  with  normal  distribution  and  standard  deviation  of  0.17,  to 
the  simulated  data. 

Figure  8.4  is  a  picture  of  the  result  of  the  reconstruction  of  the  simulated  data 
from  10  sets  of  9  frequencies  equally  spaced  in  the  range  from  1  to  4  *H. 
Random  noise  with  normal  distribution  and  standard  deviation  of  0.1%  was  again  added 
to  the  data  as  for  the  data  of  Figure  8.3,  but  here  the  first  norm  procedure  of  Eqs. 

(6  2)  was  used  on  the  predictions  to  obtain  more  accurate  predictions. 
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FIGURE  8.1:  Digitized  Image  of  Test  Object 

The  test  object  consisted  of  four  Gaussian  distributions  as  shown  within 
the  lower  box.  The  brightened  line  through  the  center  of  the  objects  indicate 
the  location  from  which  the  profile  (plotted  in  the  upper  box)  was  taken.  The 
separations  of  the  Gaussian  objects  are  3.75.  2.5.  and  1.875  inn.  and  the  full 
width  at  half  maximum  is  1.0mm  for  each  of  the  objects. 


FIGURE  8.2:  Reconstruction  of  Test  Object  from 
Rational  Function  Prediction  in  the 
Absence  of  Noise. 

The  reconstructed  object  Is  almost  of  the  same  quality  as  the  original 
test  object.  Scattering  data  were  generated  by  application  of  the  Rytov 
approximation  scattering  integral.  Projection  data  were  made  from  the  rational 
function  prediction  using  nine  equally  spaced  frequencies  from  1  to  4  MHz.  Re¬ 
construction  from  the  extrapolated  data  was  done  using  the  ART  algorithm  along 
straight  lines. 


FIGURE  8.3:  Reconstruction  of  Test  Object  in 

the  Presence  of  Random  Noise  Without 
Data  Condition. 

Random  noise  with  a  normal  distribution  of  mean  zero  and  a  standard 
deviation  of  0.1%  was  added  to  the  simulated  data  at  each  frequency.  The 
rational  function  extrapolation  was  then  obtained  from  this  noisy  data  and, 
the  reconstruction  was  carried  out  in  the  same  manner  as  for  Fig.  8.2. 


FIGURE  8.4:  Reconstruction  of  Test  Object  in  the 
Presence  of  Noise  Using  the  Data  Con¬ 
ditioning  Algorithm  of  Equation  (6.2). 

Scattering  data  was  simulated  for  10  different  sets  of  9  uniformly 
spaced  frequencies  from  1-4  MHx  and  random  noise  with  a  mean  of  zero  and  a 
standard  deviation  of  0.12  was  added  to  all  of  the  data.  The  rational 
function  extrapolation  was  applied  to  each  of  the  ten  sets  of  frequencies, 
and  then  the  first  norm  procedure  of  Eq(6.2)  was  applied  to  the  extrapola¬ 
tions.  Finally,  the  reconstruction  algorithm  was  applied  to  the  averaged 
data  in  the  same  manner  as  the  previous  two  figures. 
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for  ray  refraction.  Finally,  we  describe  an  algorithm  for  circumventing  the 
effects  of  noisy  data.  Test  illustrations  are  given,  an  simulations  of  the 
Rytov  model  corresponding  to  (*) ,  namely,  for 


-  a  **  1/  *  Jl  f  *  • 


where  L  denotes  the  straight  line  path  from  r#  to  r^  ;  we  illustrate 
the  reconstruction  of  f  in  this  case.  . 
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